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Abstract 

We study the Emptiness Formation Probability (EFP) for the spin 1/2 XXZ spin 
chain. EFP P(n) detects a formation of ferromagnetic string of the length n in the ground 
state. It is expected that EFP decays in a Gaussian way for large strings P{n) ~ n~ 7 C~ n . 
Here, we propose the explicit expressions for the rate of Gaussian decay C as well as for 
the exponent 7. In order to confirm the validity of our formulas, we employed an ab initio 
simulation technique of the density-matrix renormalization group to simulate XXZ spin 
chain of sufficient length. Furthermore, we performed Monte-Carlo integration of the 
Jimbo-Miwa multiple integral for P(n). Those numerical results for P(n) support our 
formulas fairly definitely. 



We consider interacting spin 1/2 on the one-dimensional infinite lattice. The model is 
described by the Hamiltonian 



{• s 'j' s T- 1 + s ] s j+i + A ( s j s Ui ~ 4) } ' ^ 



where S = a/2 and a are Pauli matrices. It exhibits diverse physics with varying the 
anisotropy parameter A. For A < — 1 the ground state has a long-range ferromagnetic 
order, and a finite excitation gap opens in the low-lying spectrum. For A > 1 the ground 
state develops antiferromagnetic order, the gap also exists. For moderate anisotropy 
— 1 < A < 1, the antiferromagnetic order is dissolved by quantum fluctuations. The gap 
closes, and the long-range asymptotic of correlation function (SiSj) decays as a power law 
at zero temperature. 

In this article we study the Emptiness Formation Probability (EFP), 

P(n) = (11(5+2)}' ( 2 ) 

in the domain of critical phase. The importance of EFP was emphasized in the book ||. 
It reflects the nature of the ground-state and it is a good indicator of the phase-separation. 

In spite of its seemingly complicated definition, EFP has emerged quite naturally in 
the course of studies of the quantum integrability of the XXZ spin chain. The first 
progress in calculating of EFP was achieved in Ref. [|J, where P(3) was obtained for 
the isotropic antiferromagnetic chain (A = 1). Subsequently, Jimbo and Miwa derived 
multiple integral representation for all correlation functions by means of vertex operators 
approach fi |«|. Their integral formula for EFP is: 

p(n) = {-u)-^ r^i... r^n sinh(Aa ~ Afe) 



2tt J -00 27r sinh ((A a — A& — ivr) v) 



a>b 

A sinh"- fc ((X k + ivr/2) v) sinh^ 1 ((A fc - jg/2) v) 

where the parameter v is related to the anisotropy A as 

v= icos^fA). (4) 

7T 

The integral representation was reproduced in the framework of algebraic Bethe ansatz 
in 0]. A good deal of new developments have been reported for EFP. For example, in the 
case of A = 1/2 the simple formula for P{n) has been conjectured in [JO]] and then proved 



in [11]. For A = the asymptotic form of P(n) as n > 1 was found in |5[. In either case 



the large-n asymptotic behavior is given by P(n) ~ n ^ C n with 

C = V2, 7 = i, (^ = 1/2), 

Moreover, P(n) was calculated recently for A = 1 for strings of the length n < 6, see 
papers @, §, Their result suggests that P(n) at A = 1 decays in a Gaussian way 
as well. Meanwhile, the asymptotic behavior of EFP for the whole critical regime was 
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nalyzed in field-theoretical framework, by Abanov and Korepin [12|. The Gaussian decay 
naturally appears in this approach. 

Based on those developments, it would appear reasonable that the following asymptotic 
form of EFP 

P(n)~An-^C- n \ (6) 

holds for all over the critical regime (—1 < A < 1). We propose the explicit expressions 
for the rate of Gaussian decay C and the power-law exponent 7: 

_ r 2 (l/4) J r dt sinh 2 Q)e- f 1 
( vrv 7 ^ eXP \~io Tcosh(2tz^)sinh(t) J ' ^ ^ 

1 ^ 2 

7 = 12 + 3(1^)- (8) 

One can confirm that above formulas reproduce the exact results for A = and A = 1/2 
©■ 

In order to confirm the validity of the formulas for general anisotropies < v < 1, we 
have performed extensive numerical calculations of two kinds: One is the first-principle 



simulation method, namely, the density-matrix renormalization group (DMRG) [13, 14 
|l5|| , and the other is the Monte-Carlo numerical integration of the multiple integration 
formula (||) of Jimbo and Miwa. As noted afterwards, those methods are compensative, 
and we were able to perform reliable simulations for over all the critical regime. As 
a result, we could confirm fairly definitely that the above general formulas are indeed 
correct. 

Let us turn to addressing the numerical-simulation results. In Tables 1-4, we have 
listed the DMRG results of EFP for v = 0.2, ...,0.8. (The methodological details will be 
explained afterwards.) 

In addition, in each table, we have presented the logarithm of the ratio of two adjacent 
EFPs, which should behave in the form 

ln f - 7 Ml + -) + (2n + 1) In C, (9) 

\P(n + I) J n 

according to the long-distance asymptotic formula (||). Note that, after taking the ratio, 
we are able to kill the contribution of the constant factor A of (|6|) , for which, at present, 
we have no analytical prediction. (We, however, could estimate A numerically from our 
DMRG data as is shown in Table 5.) In this way, the resultant processed data can be 
directly comparable with the analytical conjecture. 

We see that our processed DMRG data ]n(P(n)/P(n + 1)) are extremely close to the 
analytical prediction @ . Actually we find that they coincide up to two or three digits for 
general values of n and v. Therefore we conclude that our formulas for the asymptotic 
form are valid indeed for all over the regime < v < 1. Especially see the Table 6, where 
we compare the asymptotic formula with the known exact values for A = 1 |7|, ||, [| • Note 
that, in this case, we have C = r 2 (l/4)/(vr v / 27f) = 1.66925..., 7 = 1/12, and A = 0.841. 

We, however, remark that rather large discrepancies are seen for such cases either very 
large n for small v or small n for large v. The former deviations are merely due to the 
numerical round off errors. Since in computers, real numbers are stored in 8-byte space, 
and the precision is of the order of 10~ 15 at best. Hence it is in principle difficult to 
calculate the correlations less than the magnitude < 10~ 12 reliably. The latter deviations 
are not so surprising, because our general formulas should be justified for long distances of 
EFP. This short-range deviation will be further exploited in the succeeding Monte-Carlo 
numerical integration analyses. 

So far, we had analyzed processed EFP data. In Fig. 1, we present the raw (unpro- 
cessed) EFP data obtained with the DMRG method. The dotted lines are the analytical 
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Table 1: 


DMRG data for P{n) 


: v = 0.2 


n 


P(n) 


]n(P(n)/P(n + 1)) 


Asymptotics (|9|) 


1 


0.5 


1.505 


1.512 


2 


1.111 x 10" 1 


2.447 


2.447 


3 


9.614 x 10~ 3 


3.397 


3.397 


4 


3.218 x 10~ 4 


4.353 


4.353 


5 


4.140 x 10 D 


5.312 


5.312 


6 


2.041 x lO" 8 


6.276 


6.271 


7 


3.837 x 10- 11 


7.187 


7.232 


8 


2.901 x 10- 14 


4.683 


8.192 


9 


2.682 x 10~ 16 


3.923 


9.155 




Table 2: 


DMRG data for P(n) 


: v = 0.4 


n 


P(n) 


\n(P(n)/P(n + 1)) 


Asymptotics (Q) 


1 


0.5 


1.319 


1.320 


2 


1.337 x 10" 1 


2.070 


2.071 


3 


1.687 x 10~ 2 


2.851 


2.851 


4 


9.752 x 10~ 4 


3.640 


3.640 


5 


2.561 x 10 


A A O A 

4.434 


4.433 


6 


3.039 x 10~ 7 


5.231 


5.230 


7 


1.626 x 10~ 9 


6.026 


6.025 


8 


3.924 x 10" 12 


5.862 


6.823 


9 


1.117 x 10~ 14 


2.824 


7.621 




Table 3: 


DMRG data for P(n) 


: v = 0.6 


n 


P(n) 


\a(P(n)/P(n + 1)) 


Asymptotics @ 


1 


0.5 


1.106 


1.125 


2 


1.655 x lO^ 1 


1.594 


1.588 


3 


3.360 x 10~ 2 


2.121 


2.115 


4 


4.028 x 10~ 3 


2.667 


2.663 


5 


2.798 x 10 


6.226 


3.221 


6 


1.115 x 10~ 5 


3.784 


3.783 


7 


2.534 x 10~ 7 


4.348 


4.348 


8 


3.278 x 10~ 9 


4.909 


4.915 


9 


2.420 x 10" 11 


5.231 


5.483 
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Table 4: DMRG data for P(n) 


: u = 0.8 


n 


P{n) 


\n(P(n)/P{n + 1)) 


Asymptotics (0) 


1 


0.5 


0.896 


1.250 


2 


2.041 x 10" 1 


1.120 


1.221 


3 


6.663 x 10" 2 


1.361 


1.387 


4 


1.709 x lO- 2 


1.615 


1.615 


5 


3.400 x 10~ 3 


1.879 


1.870 


6 


5.196 x 10~ 4 


2.150 


2.139 


7 


6.050 x 10~ 5 


2.430 


2.417 


8 


5.336 x 10~ 6 


2.708 


2.701 


9 


3.549 x HT 7 


2.996 


2.988 



Table 5: Numerical estimation of the factor A 



V 


A 


C 


7 


A 


0.0 


1.0 


1.66925 


0.0833 


0.841 


0.2 


0.8090 


1.61803 


0.1 


0.816 


0.4 


0.3090 


1.49207 


0.1722 


0.747 


0.6 


-0.3090 


1.33168 


0.3833 


0.68 


0.8 


-0.8090 


1.16287 


1.15 


0.9 



Table 6: Comparison with exact values for P(n) at A = 1 



n 


P(n) 




\n(P(n)/P(n + 1)) 


Asymptotics (H) 


1 


0.5 




1.58685 


1.59489 


2 


1.02284 x 10" 


-i 


2.59643 


2.59566 


3 


7.62415 x 10" 


-3 


3.60989 


3.61060 


4 


2.06270 x 10- 


-4 


4.63019 


4.62998 


5 


2.01172 x 10- 


-6 


5.65115 


5.65133 


6 


7.06812 x 10' 


-9 
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conjecture with the factors given in Table 5. We see that our first-principle data are 
well-fitted by our general formulas (^) and (||) for various values of v. However, as noted 
above, for some cases, there appear deviations. 




Figure 1: EFP is plotted for various valued of distances n and anisotropies v. We have employed 
the DMRG method. The dotted lines are our conjecture @ with the formulas (0) and (H). 



Let us turn to addressing details of the DMRG method [13, 14, 15]. DMRG is a 
sort of computer-aided real-space decimation, where the number of states for a renormal- 
ized block is retained within a bound tractable in computers. In this way, through the 
successive applications of DMRG, we are able to access very large system sizes. Hence, 
this method is suitable for surveying the long-distance asymptotic behavior of EFP. As a 
matter of fact, efficiency of the method was demonstrated in the preceding works ||: 
This success was rather unexpected, because in general, DMRG is not very efficient for 
such systems exhibiting criticality. In particular, long-distance asymptotic behavior of 
two-point correlation function is deteriorated rather severely. Those pathologies would 
be due to the truncation of bases through numerical renormalizations. However, as far 
as EFP is concerned, DMRG is proven to be free from such difficulties. Below, we will 
overview technical points of our numerical simulation. Our algorithm is standard [14|] 



and we refer readers to a recent proceeding [15] for full account of methodological details. 



We have employed the so-called infinite algorithm, which is adequate to investigate bulk 
properties at the ground state. (For the purpose of studying finite-size scaling behavior, 
the finite algorithm would be more suitable.) We have remained, at most, m = 300 bases 
for a renormalized block. The density matrix eigenvalue {w a } of remaining bases indi- 
cates the statistical weight of each remained state. We found w a > 5 • 10~ 12 : That is, we 
have remained almost all relevant states with appreciable statistical weight w a > 5- 10~ 12 , 
which may indicate error of the present simulation. We have repeated 300 renormaliza- 
tions, and hence, total system size extends to L = 600. The DMRG data alternates in 
turn through renormalizations: Note that the number of spins consisting a renormalized 
block increases by one after another through renormalizations. The problem is that the 
Hilbert-space structures are incompatible with respect to those cases whether the block 
contains even or odd spins. (Haldane system (S = 1) is not affected by this difficulty.) 
Therefore, we have taken arithmetic mean over those two cases. 

The DMRG method is not efficient in the close vicinity of the ferromagnetic isotropic 
point v = 1. Because in the vicinity of that point, the spin-wave velocity tends to 
vanish, and there appear numerous nearly-degenerated low-lying levels. Those nearly- 
degenerated levels are very hard to resolve in the process of numerical diagonalization. 
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The diagonalization is a significant part in the DMRG procedure, and hence, DMRG 
becomes hardly applicable. In order to compensate this drawback, we have employed 
another numerical method, that is, the Monte-Carlo integration. We will explain it in the 
following. 

r=— ^ 1 1 1 1 1 



-12 



= 0.80 
0.90 
0.95 
0.99 

(1/2)" 
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Figure 2: EFP is plotted for various valued of distances n and anisotropies v. We have employed 
the Monte-Carlo integration method for Eq. (^). The dotted line is Eq. (^) for v = 0.8. 

In Fig. 2, we have presented numerical data with Monte-Carlo integration for the 
multiple integral formula (||) of Jimbo and Miwa. The technical details will be explained 
afterwards. From the plot, we see that the decay of EFP is gradually modified as for 
!/—►]., and surprisingly enough, the decay approaches the simple exponential formula 
of (l/2) n eventually. This feature is precisely in accord with the aforementioned finding 
that DMRG data deviate from the Gaussian formula as for v — > 1. Therefore, we see that 
at least for short-distance asymptotic form of EFP is governed by the simple exponential 
decay 0.5™ as u — ► 1. As a matter of fact, right at v = 1, one may easily verify the pure 
exponential decay P(n) from the Jimbo-Miwa formula. 

The crossover from Gaussian to exponential decay suggests that the quantum fluctu- 
ations are suppressed as v — * 1. As as matter of fact, field theoretical consideration 
reveled that EFP measures the probability of formations of ferromagnetic islands with size 
n surrounded by antiferromagnetic background from the viewpoint of Euclidean space- 
time. Therefore, our result suggests that this Euclidean space-time picture is deteriorated 
for v — > 1, because the spin- wave velocity tends to vanish, and thus, the size of the islands 
grows along the imaginary-time direction abruptly. 

In the following, we explain technical details of the Monte-Carlo integration. We 



have used the subroutine package described in the textbook [16]. As a random number 
generator, we have employed "Mersenne Twistor", MT19937 in Ref. [jT^ ]. For each plot, 
we have performed eight-billion Monte-Carlo steps. This main Monte-Carlo procedure is 
preceded by preliminary five million Monte-Carlo steps which are aimed to improve the 
efficiency of Monte-Carlo sampling by surveying the integrand in the multi-dimensional 
space. Irrespective of n and u, we found that the statistical errors are of the order of 

nr 8 . 

To summarize, we have advocated compact explicit formulas (0) and (H) for the long- 
distance asymptotic form behavior of EFP. The formulas reproduce the presently-available 
exact results obtained at special solvable points of A = and 1/2, and therefore, we 
expect that our formulas are valid for over all the critical regime. We have performed 
extensive simulations of DMRG and Monte-Carlo integration. As a result, we found that 
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our general formulas are indeed correct over the critical regime; see Tables and Fig. 1. In 
other words, EFP decays obeying the Gaussian formula (||) in the domain of criticality, 
and the decay rate becomes slower in the ferromagnetic side. This tendency is understood 
by a physical picture that EFP detects a formation of condensed particles (ferromagnetic 
string). Apparently, for ferromagnetic side, the condensation becomes promoted, and EFP 
gets enhanced, albeit the asymptotic form is still maintained to be Gaussian. However, as 
for the extreme limit to the ferromagnetic isotropic point, the Gaussian decay is smeared 
out by an exponential form at least for short distances. This fact reflects the suppression 
of quantum fluctuations in the ferromagnetic side. In other words, this crossover may 
be regarded as the precursor of the onset of true (long-range order) phase-separation for 
A < — 1, where EFP should not decay, but may saturate to a certain finite resident value. 
Finally we like to remark that the rate of Gaussian decay C, Eq.(^), can be evaluated 
analytically for some other rational values of v. For example, we have found 

1 G 

C\ v= \u = 24 e~, (G : Catalan number), 

o 1^=1/6 = iy-, c 1^=2/3- — ^= — , (10) 

by use of the integral representations for the logarithm of (multiple) Gamma functions. 
These facts may indicate a possibility to prove our analytic formulas (^) rigorously for 
these values of v just in the similar way as v = 1/2 and 1/3. 
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